!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>


!     subroutine LIMBES_boundary computes boundary condition which a combination
!     of diffuse and reflexive boudary condition (i.e. Maxwell's)
!     used common variables:
!     SIMUL_wall_x(1:LIMBES_nmax),SIMUL_wall_y(1:LIMBES_nmax)
!     -> node positions xi,yi (integer*8)
!      SIMUL_wall_ex(1:LIMBES_nmax),SIMUL_wall_ey(1:LIMBES_nmax)
!     -> normal vector component (ex,ey) to the wall at the node (integer*8)
!     SIMUL_wall_alpha(1:LIMBES_nmax) -> thermal accomodation coefficient alpha (real*8)
!     SIMUL_wall_u(1:LIMBES_nmax) -> wall tangentiel velocity uw at this node (real*8)
!     SIMUL_wall_RT(1:LIMBES_nmax) -> wall temperature RTw at this node (real*8)
!     all these wall information are read in the file ../data/boundary.in with this format:
!     SIMUL_wall_nodesnumber
!     ni xi yi ex ey alpha uw RTw
!     needed argument:
!     integer*8 ni node index

subroutine LIMBES_boundary_1wall(nwi)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: ni,nwi,vi,nv,RLUD
  real*8 :: vx,vy,alpha,uw,ux,uy,RTw,ex,ey,rho_w,P_w,Qw
  real*8 :: F_right(1:LIMBES_nv)
  real*8 :: F_left(1:LIMBES_nv)
  real*8 :: F_up(1:LIMBES_nv)
  real*8 :: F_down(1:LIMBES_nv) 
  
  nv=LIMBES_nv
  ni=SIMUL_1wall_cell(nwi)
 
  !    Fluxes computation on boundaries cells 
  RLUD=1
  call F_flux(ni,F_right,RLUD)
  RLUD=2
  call F_flux(ni,F_left,RLUD)
  RLUD=3
  call F_flux(ni,F_up,RLUD)
  RLUD=4
  call F_flux(ni,F_down,RLUD)

  select case (SIMUL_1wall_boundary_type(nwi))
  case("S") ! S for Specular boundary condition       
     !     Fluxes computation on boundaries cells
     do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then
           if (RLUD.EQ.1) then
              call Normal_vector(ex,ey,RLUD)
              call F_specular_flux(F_right,ex,ey)                 
           endif
           if (RLUD.EQ.2) then
              call Normal_vector(ex,ey,RLUD)
              call F_specular_flux(F_left,ex,ey)                
           endif
           if (RLUD.EQ.3) then
              call Normal_vector(ex,ey,RLUD)
              call F_specular_flux(F_up,ex,ey)                
           endif
           if (RLUD.EQ.4) then
              call Normal_vector(ex,ey,RLUD)
              call F_specular_flux(F_down,ex,ey)                
           endif
        endif
     enddo
     !     --------------------------------------------------------------------          
  case("T") ! T for Thermal boundary condition
     alpha=SIMUL_1wall_alpha(nwi,1)
     uw=SIMUL_1wall_alpha(nwi,2)
     RTw=SIMUL_1wall_alpha(nwi,3)*FLUID_R
     !   Fluxes computation on boundaries cells
     do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then
           if (RLUD.EQ.1) then
              call Normal_vector(ex,ey,RLUD)
              call diffuse_boundary(alpha,uw,RTw,ex,ey,F_right)
           endif
           if (RLUD.EQ.2) then
              call Normal_vector(ex,ey,RLUD)
              call diffuse_boundary(alpha,uw,RTw,ex,ey,F_left)
           endif
           if (RLUD.EQ.3) then
              call Normal_vector(ex,ey,RLUD)
              call diffuse_boundary(alpha,uw,RTw,ex,ey,F_up)
           endif
           if (RLUD.EQ.4) then
              call Normal_vector(ex,ey,RLUD)
              call diffuse_boundary(alpha,uw,RTw,ex,ey,F_down)
           endif
        endif
     enddo
     !     --------------------------------------------------------------------          
  case("H") ! H for Heat flux boundary condition
     alpha=SIMUL_1wall_alpha(nwi,1)
     uw=SIMUL_1wall_alpha(nwi,2)
     Qw=SIMUL_1wall_alpha(nwi,4)
     !   Fluxes computation on boundaries cells
     do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then
           if (RLUD.EQ.1) then
              call Normal_vector(ex,ey,RLUD)
              call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_right)
           endif
           if (RLUD.EQ.2) then
              call Normal_vector(ex,ey,RLUD)
              call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_left)
           endif
           if (RLUD.EQ.3) then
              call Normal_vector(ex,ey,RLUD)
              call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_up)
           endif
           if (RLUD.EQ.4) then
              call Normal_vector(ex,ey,RLUD)
              call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_down)
           endif
        endif
     enddo
     SIMUL_1wall_alpha(nwi,3)=RTw/FLUID_R
     !     --------------------------------------------------------------------         
  case("P") ! P for Pressure boundary conditions
     P_w=SIMUL_1wall_alpha(nwi,1)
     do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then
           if (RLUD.EQ.1) then
              call pressure_boundary(ni,P_w,RLUD,F_right)
           endif
           if (RLUD.EQ.2) then
              call pressure_boundary(ni,P_w,RLUD,F_left)
           endif
           if (RLUD.EQ.3) then
              call pressure_boundary(ni,P_w,RLUD,F_up)
           endif
           if (RLUD.EQ.4) then
              call pressure_boundary(ni,P_w,RLUD,F_down)
           endif
        endif
     enddo
     !     --------------------------------------------------------------------         
  case("PT") ! PT for Pressure boundary conditions and imposed Temperature
     P_w=SIMUL_1wall_alpha(nwi,1)
     RTw=SIMUL_1wall_alpha(nwi,2)*FLUID_R
     do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then
           if (RLUD.EQ.1) then
              call pressureT_boundary(ni,P_w,RTw,RLUD,F_right)
           endif
           if (RLUD.EQ.2) then
              call pressureT_boundary(ni,P_w,RTw,RLUD,F_left)
           endif
           if (RLUD.EQ.3) then
              call pressureT_boundary(ni,P_w,RTw,RLUD,F_up)
           endif
           if (RLUD.EQ.4) then
              call pressureT_boundary(ni,P_w,RTw,RLUD,F_down)
           endif
        endif
     enddo
     !     --------------------------------------------------------------------         
  case("F") ! F for Flow boundary conditions
     rho_w=SIMUL_1wall_alpha(nwi,1)
     ux=SIMUL_1wall_alpha(nwi,2)
     uy=SIMUL_1wall_alpha(nwi,3)
     RTw=SIMUL_1wall_alpha(nwi,4)*FLUID_R
     do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then
           if (RLUD.EQ.1) then
              call Normal_vector(ex,ey,RLUD)
              call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_right)
           endif
           if (RLUD.EQ.2) then
              call Normal_vector(ex,ey,RLUD)
              call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_left)
           endif
           if (RLUD.EQ.3) then
              call Normal_vector(ex,ey,RLUD)
              call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_up)
           endif
           if (RLUD.EQ.4) then
              call Normal_vector(ex,ey,RLUD)
              call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_down)
           endif
        endif
     enddo
     !     --------------------------------------------------------------------         
  case default
     print*,'error from LIMBES_boundary_1wall boundary condition type not known'
     print*,'SIMUL_1wall_boundary_type(ni)=', SIMUL_1wall_boundary_type(nwi),nwi,ni
     stop
  end select
  
  !     v.grad(f)
  LIMBES_vgradf(ni,:)=LIMBES_vx(:)*(F_right(:) - F_left(:))/SIMUL_dx + LIMBES_vy(:)*(F_up(:) - F_down(:))/SIMUL_dx
  !     f at the wall
  do RLUD=1,4
        if  (SIMUL_1wall_RLUD(nwi,RLUD).EQ.1) then 
           if (RLUD.EQ.1) LIMBES_f_1wall(nwi,:)=F_right(:)
           if (RLUD.EQ.2) LIMBES_f_1wall(nwi,:)=F_left(:)
           if (RLUD.EQ.3) LIMBES_f_1wall(nwi,:)=F_up(:)
           if (RLUD.EQ.4) LIMBES_f_1wall(nwi,:)=F_down(:)
        endif
     enddo
     

end subroutine LIMBES_boundary_1wall
      
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine LIMBES_boundary_2wall(nwi)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: ni,nwi,vi,nv,RLUD,nwii,wall_RLUD
  real*8 :: alpha,ux,uy,uw,RTw,rho_w,P_w,Qw
  real*8 :: vx,vy,ex,ey
  real*8 :: F_right(1:LIMBES_nv)
  real*8 :: F_left(1:LIMBES_nv)
  real*8 :: F_up(1:LIMBES_nv)
  real*8 :: F_down(1:LIMBES_nv) 
  logical :: test

  nv=LIMBES_nv
  ni=SIMUL_2wall_cell(nwi)  

  !     Fluxes computation on boundaries cells 
  RLUD=1
  call F_flux(ni,F_right,RLUD)
  RLUD=2
  call F_flux(ni,F_left,RLUD)
  RLUD=3
  call F_flux(ni,F_up,RLUD)
  RLUD=4
  call F_flux(ni,F_down,RLUD)



  wall_RLUD=0
  do nwii=1,2
     test=.true.
     select case (SIMUL_2wall_boundary_type(nwi,nwii))
     case("S") ! S for Specular boundary condition         
        !     Diffuse fluxes on boundary
        do RLUD=1,4
           if (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then
              if (RLUD.EQ.1.AND.wall_RLUD.NE.1.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call F_specular_flux(F_right,ex,ey)
              endif
              if (RLUD.EQ.2.AND.wall_RLUD.NE.2.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call F_specular_flux(F_left,ex,ey)
              endif
              if (RLUD.EQ.3.AND.wall_RLUD.NE.3.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call F_specular_flux(F_up,ex,ey)
              endif
              if (RLUD.EQ.4.AND.wall_RLUD.NE.4.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call F_specular_flux(F_down,ex,ey)
              endif
           endif
        enddo
        !     --------------------------------------------------------------------   
     case("T") ! T for Thermal boundary condition   
        if (nwii.EQ.1) then
           alpha=SIMUL_2wall_alpha1(nwi,1)
           uw=SIMUL_2wall_alpha1(nwi,2)
           RTw=SIMUL_2wall_alpha1(nwi,3)*FLUID_R
        endif
        if (nwii.EQ.2) then
           alpha=SIMUL_2wall_alpha2(nwi,1)
           uw=SIMUL_2wall_alpha2(nwi,2)
           RTw=SIMUL_2wall_alpha2(nwi,3)*FLUID_R
        endif
        do RLUD=1,4
           if (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then     
              if (RLUD.EQ.1.AND.wall_RLUD.NE.1.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call diffuse_boundary(alpha,uw,RTw,ex,ey,F_right)
              endif
              if (RLUD.EQ.2.AND.wall_RLUD.NE.2.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call diffuse_boundary(alpha,uw,RTw,ex,ey,F_left)
              endif
              if (RLUD.EQ.3.AND.wall_RLUD.NE.3.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call diffuse_boundary(alpha,uw,RTw,ex,ey,F_up)
              endif
              if (RLUD.EQ.4.AND.wall_RLUD.NE.4.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call diffuse_boundary(alpha,uw,RTw,ex,ey,F_down)
              endif
           endif
        enddo

	

        !     --------------------------------------------------------------------   
     case("H") ! H for Heat flux boundary condition  
        if (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then
           if (nwii.EQ.1) then
              alpha=SIMUL_2wall_alpha1(nwi,1)
              uw=SIMUL_2wall_alpha1(nwi,2)
              Qw=SIMUL_2wall_alpha1(nwi,4)
           endif
           if (nwii.EQ.2) then
              alpha=SIMUL_2wall_alpha2(nwi,1)
              uw=SIMUL_2wall_alpha2(nwi,2)
              Qw=SIMUL_2wall_alpha2(nwi,4)
           endif
           do RLUD=1,4
              if (RLUD.EQ.1.AND.wall_RLUD.NE.1.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_right)
              endif
              if (RLUD.EQ.2.AND.wall_RLUD.NE.2.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_left)
              endif
              if (RLUD.EQ.3.AND.wall_RLUD.NE.3.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_up)
              endif
              if (RLUD.EQ.4.AND.wall_RLUD.NE.4.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_down)
              endif
           enddo
           if (nwii.EQ.1) then
              SIMUL_2wall_alpha1(nwi,3)=RTw/FLUID_R
           endif
           if (nwii.EQ.2) then
              SIMUL_2wall_alpha2(nwi,3)=RTw/FLUID_R
           endif
        endif
        !     --------------------------------------------------------------------
     case("P")  ! P for Pressure boundary conditions  
        if (nwii.EQ.1) then
           P_w=SIMUL_2wall_alpha1(nwi,1)     
        endif
        if (nwii.EQ.2) then
           P_w=SIMUL_2wall_alpha2(nwi,1)
        endif
        do RLUD=1,4
           if  (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then
              if (RLUD.EQ.1.AND.wall_RLUD.NE.1.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressure_boundary(ni,P_w,RLUD,F_right)
              endif
              if (RLUD.EQ.2.AND.wall_RLUD.NE.2.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressure_boundary(ni,P_w,RLUD,F_left)
              endif
              if (RLUD.EQ.3.AND.wall_RLUD.NE.3.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressure_boundary(ni,P_w,RLUD,F_up)
              endif
              if (RLUD.EQ.4.AND.wall_RLUD.NE.4.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressure_boundary(ni,P_w,RLUD,F_down)
              endif
           endif
        enddo
        !     --------------------------------------------------------------------
     case("PT")  ! P for Pressure boundary conditions and imposed Temperature  
        if (nwii.EQ.1) then
           P_w=SIMUL_2wall_alpha1(nwi,1)
           RTw=SIMUL_2wall_alpha1(nwi,2)*FLUID_R
        endif
        if (nwii.EQ.2) then
           P_w=SIMUL_2wall_alpha2(nwi,1)
           RTw=SIMUL_2wall_alpha2(nwi,2)*FLUID_R
        endif
        do RLUD=1,4
           if  (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then
              if (RLUD.EQ.1.AND.wall_RLUD.NE.1.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressureT_boundary(ni,P_w,RTw,RLUD,F_right)
              endif
              if (RLUD.EQ.2.AND.wall_RLUD.NE.2.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressureT_boundary(ni,P_w,RTw,RLUD,F_left)
              endif
              if (RLUD.EQ.3.AND.wall_RLUD.NE.3.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressureT_boundary(ni,P_w,RTw,RLUD,F_up)
              endif
              if (RLUD.EQ.4.AND.wall_RLUD.NE.4.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call pressureT_boundary(ni,P_w,RTw,RLUD,F_down)
              endif
           endif
        enddo
        !     --------------------------------------------------------------------
     case("F")  ! F for Flow boundary conditions  
        if (nwii.EQ.1) then
           rho_w=SIMUL_2wall_alpha1(nwi,1)
           ux=SIMUL_2wall_alpha1(nwi,2)
           uy=SIMUL_2wall_alpha1(nwi,3)
           RTw=SIMUL_2wall_alpha1(nwi,4)*FLUID_R
        endif
        if (nwii.EQ.2) then
           rho_w=SIMUL_2wall_alpha2(nwi,1)
           ux=SIMUL_2wall_alpha2(nwi,2)
           uy=SIMUL_2wall_alpha2(nwi,3)
           RTw=SIMUL_2wall_alpha2(nwi,4)*FLUID_R
        endif
        do RLUD=1,4
           if  (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then
              if (RLUD.EQ.1.AND.wall_RLUD.NE.1.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_right)
              endif
              if (RLUD.EQ.2.AND.wall_RLUD.NE.2.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_left)
              endif
              if (RLUD.EQ.3.AND.wall_RLUD.NE.3.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_up)
              endif
              if (RLUD.EQ.4.AND.wall_RLUD.NE.4.AND.test) then
                 wall_RLUD=RLUD
                 test=.false.
                 call Normal_vector(ex,ey,RLUD)
                 call flow_boundary(rho_w,ux,uy,RTw,ex,ey,F_down)
              endif
           endif
        enddo
        !     --------------------------------------------------------------------
     case default
        print*,'error from LIMBES_boundary_1wall:boundary condition type not known'
        print*,'SIMUL_1wall_boundary_type(ni,nwii)=', SIMUL_2wall_boundary_type(ni,nwii)
        stop
     end select
  enddo


  !     v.grad(f)
  LIMBES_vgradf(ni,:)= LIMBES_vx(:)*(F_right(:) - F_left(:))/SIMUL_dx + LIMBES_vy(:)*(F_up(:) - F_down(:))/SIMUL_dx
  !     f at the wall
  do RLUD=1,4
     if  (SIMUL_2wall_RLUD(nwi,RLUD).EQ.1) then 
        if (RLUD.EQ.1) LIMBES_f_2wall(nwi,:,1)=F_right(:)
        if (RLUD.EQ.2) LIMBES_f_2wall(nwi,:,1)=F_left(:)
        if (RLUD.EQ.3) LIMBES_f_2wall(nwi,:,2)=F_up(:)
        if (RLUD.EQ.4) LIMBES_f_2wall(nwi,:,2)=F_down(:)
     endif
  enddo
  
end subroutine LIMBES_boundary_2wall
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine F_specular_flux(F_RLUD,ex,ey)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: vi,vi1,nv
  real*8 :: ex,ey,vdote,vx,vy
  real*8 :: F_RLUD(1:LIMBES_nv)
     
  nv=LIMBES_nv
      
  !     specular boundary condition
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then           
        if (dabs(ex).GT.0) vi1=LIMBES_opposite_vx(vi)
        if (dabs(ey).GT.0) vi1=LIMBES_opposite_vy(vi)
        F_RLUD(vi)=F_RLUD(vi1)
     endif
  enddo
  
end subroutine F_specular_flux

!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine diffuse_boundary(alpha,uw,RTw,ex,ey,F_bound)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: vi,nv
  real*8 :: alpha,uw,RTw,ux,uy
  real*8 :: ex,ey,vdote,vx,vy,vn,w_f,w_feq,J_in,J_out,rho_w
  real*8 :: F_bound(1:LIMBES_nv) 
  real*8 :: feq(1:LIMBES_nv) 
  
  nv=LIMBES_nv
  
  !    specular reflexion
  call F_specular_flux(F_bound,ex,ey)

  J_in=zero
  J_out=zero
  ux=dabs(ey)*uw
  uy=dabs(ex)*uw      
  
  call LIMBES_calc_feq(one,ux,uy,RTw,feq)     
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then
        vn=vx*dabs(ex) + vy*dabs(ey)  
        w_f=LIMBES_w(vi)*F_bound(vi)
        w_feq=LIMBES_w(vi)*feq(vi)                                 
        J_in=J_in + vn*w_f                               
        J_out=J_out + vn*w_feq
     endif
  enddo
  rho_w=J_in/J_out
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then 
        F_bound(vi)=alpha*rho_w*feq(vi)+(one-alpha)*F_bound(vi)
     endif
  enddo

end subroutine diffuse_boundary
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine heat_boundary(alpha,uw,RTw,Qw,ex,ey,F_bound)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: vi,nv
  real*8 :: alpha,uw,RTw,ux,uy
  real*8 :: ex,ey,vdote,vx,vy,vn,w_f,w_feq,J_in,J_out,rho_w,Qw
  real*8 :: F_bound(1:LIMBES_nv) 
  real*8 :: feq(1:LIMBES_nv) 
  
  nv=LIMBES_nv
  
  !    specular reflexion
  call F_specular_flux(F_bound,ex,ey)

  J_in=zero
  J_out=zero
  ux=dabs(ey)*uw
  uy=dabs(ex)*uw      
  
  call LIMBES_calc_feq(one,ux,uy,RTw,feq)     
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then
        vn=vx*dabs(ex) + vy*dabs(ey)  
        w_f=LIMBES_w(vi)*F_bound(vi)
        w_feq=LIMBES_w(vi)*feq(vi)                                 
        J_in=J_in + vn*w_f                               
        J_out=J_out + vn*w_feq
     endif
  enddo
  rho_w=J_in/J_out
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then 
        F_bound(vi)=alpha*rho_w*feq(vi)+(one-alpha)*F_bound(vi)
     endif
  enddo

end subroutine heat_boundary

!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine pressure_boundary(ni,P_w,RLUD,F_bound)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: vi,ni,nv,xi,yi,n2,RLUD
  real*8 :: RT,ux,uy
  real*8 :: ux1,ux2,uy1,uy2,RT1,RT2,rho1,rho2
  real*8 :: ex,ey,vdote,vx,vy,P_w,rho_w
  real*8 :: F_bound(1:LIMBES_nv),F_neq_bound(1:LIMBES_nv)
  real*8 :: F_neq_bound1(1:LIMBES_nv),F_neq_bound2(1:LIMBES_nv)
  real*8 :: feq1(1:LIMBES_nv),feq2(1:LIMBES_nv)
  real*8 :: feq(1:LIMBES_nv)
  nv=LIMBES_nv
      
  call Normal_vector(ex,ey,RLUD)
      
  xi=LIMBES_cell_xi(ni)
  yi=LIMBES_cell_yi(ni)
  if (RLUD.EQ.1) then
     n2=LIMBES_cell_ni(xi-1,yi)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  if (RLUD.EQ.2) then
     n2=LIMBES_cell_ni(xi+1,yi)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  if (RLUD.EQ.3) then
     n2=LIMBES_cell_ni(xi,yi-1)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  if (RLUD.EQ.4) then
     n2=LIMBES_cell_ni(xi,yi+1)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  
  ux=1.5d+0*ux1 - 0.5d+0*ux2
  uy=1.5d+0*uy1 - 0.5d+0*uy2
  !      RT=1.5d+0*RT1 - 0.5d+0*RT2
  RT=RT1
  rho_w=P_w/RT
  
  call LIMBES_calc_feq(rho_w,ux,uy,RT,feq)    
  
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then 
        F_neq_bound(vi)=1.5D+0*F_neq_bound1(vi) - 0.5D+0* F_neq_bound2(vi)
        F_bound(vi)=feq(vi) + F_neq_bound(vi)
     endif
  enddo
  
end subroutine pressure_boundary
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine pressureT_boundary(ni,P_w,RTw,RLUD,F_bound)
  use LIMBES_mod_fluid_config
  use LIMBES_mod_var
  implicit none
  integer*8 :: vi,ni,nv,xi,yi,n2,RLUD
  real*8 :: RTw,ux,uy
  real*8 :: ux1,ux2,uy1,uy2
  real*8 :: ex,ey,vdote,vx,vy,rho_w,P_w,rho1,rho2,RT1,RT2
  real*8 :: F_bound(1:LIMBES_nv),F_neq_bound(1:LIMBES_nv)
  real*8 :: F_neq_bound1(1:LIMBES_nv),F_neq_bound2(1:LIMBES_nv)
  real*8 :: feq1(1:LIMBES_nv),feq2(1:LIMBES_nv)
  real*8 :: feq(1:LIMBES_nv)
  
  nv=LIMBES_nv
  call Normal_vector(ex,ey,RLUD)
  
  xi=LIMBES_cell_xi(ni)
  yi=LIMBES_cell_yi(ni)
  if (RLUD.EQ.1) then
     n2=LIMBES_cell_ni(xi-1,yi)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  if (RLUD.EQ.2) then
     n2=LIMBES_cell_ni(xi+1,yi)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  if (RLUD.EQ.3) then
     n2=LIMBES_cell_ni(xi,yi-1)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif
  if (RLUD.EQ.4) then
     n2=LIMBES_cell_ni(xi,yi+1)
     ux1=LIMBES_ux(ni)
     ux2=LIMBES_ux(n2)
     uy1=LIMBES_uy(ni)
     uy2=LIMBES_uy(n2)
     rho1=LIMBES_rho(ni)
     rho2=LIMBES_rho(n2)
     RT1=LIMBES_RT(ni)
     RT2=LIMBES_RT(n2)
     call LIMBES_calc_feq(rho1,ux1,uy1,RT1,feq1)
     call LIMBES_calc_feq(rho2,ux2,uy2,RT2,feq2)
     do vi=1,nv
        F_neq_bound1(vi)=LIMBES_f(ni,vi) - feq1(vi)
        F_neq_bound2(vi)=LIMBES_f(n2,vi) - feq2(vi)
     enddo
  endif

  ux=1.5d+0*ux1 - 0.5d+0*ux2
  uy=1.5d+0*uy1 - 0.5d+0*uy2
  rho_w=P_w/RTw
  call LIMBES_calc_feq(rho_w,ux,uy,RTw,feq)   
    
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then 
        F_neq_bound(vi)=1.5D+0*F_neq_bound1(vi) - 0.5D+0* F_neq_bound2(vi)
        F_bound(vi)=feq(vi) + F_neq_bound(vi)
     endif
  enddo
  
end subroutine pressureT_boundary
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine flow_boundary(rho,ux,uy,RT,ex,ey,F_bound)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 ::  vi,nv
  real*8 :: rho,RT,ux,uy
  real*8 :: ex,ey,vdote,vx,vy
  real*8 :: F_bound(1:LIMBES_nv) 
  real*8 :: feq(1:LIMBES_nv) 
  
  nv=LIMBES_nv
  call LIMBES_calc_feq(rho,ux,uy,RT,feq)    
  
  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     vdote=vx*ex + vy*ey
     if (vdote.GT.zero) then 
        F_bound(vi)=feq(vi)
     endif
  enddo

end subroutine flow_boundary
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine Normal_vector(ex,ey,RLUD)
  use LIMBES_mod_fluid_config
  implicit none 
  real*8 :: ex,ey
  integer*8 :: RLUD

  if (RLUD.EQ.1) then
     ex=-one
     ey=zero
  endif
  if (RLUD.EQ.2) then
     ex=one
     ey=zero
  endif
  if (RLUD.EQ.3) then
     ex=zero
     ey=-one
  endif
  if (RLUD.EQ.4) then
     ex=zero
     ey=one    
  endif
  
end  subroutine Normal_vector
      
